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Abstract 

We report the results of a numerical study of nonequilibrium steady states for a class of 
Hamiltonian models. In these models of coupled matter-energy transport, particles exchange 
energy through collisions with pinned-down rotating disks. In [16], Eckmann and Young stud- 
ied ID chains and showed that certain simple formulas give excellent approximations of energy 
and particle density profiles. Keeping the basic mode of interaction in lfl6l . we extend their 
prediction scheme to a number of new settings: 2D systems on different lattices, driven by a va- 
riety of boundary (heat bath) conditions including the use of thermostats. Particle-conserving 
models of the same type are shown to behave similarly. The second half of this paper exam- 
ines memory and finite-size effects, which appear to impact only minimally the profiles of the 
models tested in lfl6l . We demonstrate that these effects can be significant or insignificant 
depending on the local geometry. Dynamical mechanisms are proposed, and in the case of di- 
rectional bias in particle trajectories due to memory, correction schemes are derived and shown 
to give accurate predictions. 

Introduction 

Explaining far-from-equilibrium macroscopic phenomena on the basis of microscopic dynamical 
mechanisms is one of the basic problems of nonequilibrium statistical mechanics. This paper 
presents the results of some numerical studies aimed at shedding light on "micro-to-macro" ques- 
tions for a class of mechanical models introduced in ll35l l27l [T6ll as paradigms for transport pro- 
cesses. Here as in lfl6l . the basic idea is that macroscopic quantities such as energy and density 
profiles can be deduced - easily and fairly accurately - from certain information on local dynamics 
provided one makes a couple of simplifying assumptions. In this paper, we extend the results in 
lfl6l to a larger class of models and at the same time examine more closely the extent to which these 
assumptions hold. These questions lead to a systematic study of how geometry impacts memory 
and finite-size effects. 
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The models studied in this paper describe a coupled transport of matter and energy. The local 
dynamics are purely deterministic and energy-conserving; that is what we mean by "Hamiltonian" 
in this paper. The system is attached to (unequal) stochastic heat baths which both inject particles 
into the system and absorb those that leave. Placed throughout the system, evenly spaced, are rotat- 
ing disks that are nailed down at their centers; particles exchange energy with them upon collisions. 
This interaction was first used in [[351 and 11271 ; it takes the place of direct particle-particle interac- 
tions, which are more delicate to control. The authors of IfTBI introduced the following conceptual 
simplification: each rotating disk is enclosed in a well-defined area called a cell, with small enough 
passageways between cells to ensure equilibration within each cell. Motivated by this "decoupling" 
of internal cell dynamics and cell-to-cell traffic, they proposed a scheme that yields simple explicit 
expressions for approximate energy profiles and particle densities. The scheme is based on two 
assumptions: that cell-to-cell movement takes essentially the form of an unbiased random walk, 
and that within each cell, the system attains local thermal equilibrium quickly. They found striking 
agreement between their proposed formulas and simulation results, with barely visible errors for 
chains consisting of no more than 30 cells. 

The aims of this paper are twofold. One is a straightforward extension of the results in lfT6l to a 
larger class of models. Keeping the basic interaction as well as cell structure in the last paragraph, 
we have found the prediction scheme in IfTBI - suitably modified - to be very flexible. To give 
a sense of the type of generalizations, our simulations are in 2D, but we expect similar results to 
hold in 3D (with rotating balls in the place of rotating disks). In dimensions greater than one, there 
are different choices of lattices and boundary conditions. We show that this scheme is equally 
valid for rectangular and hexagonal lattices, and for all of the boundary conditions tested. The 
resulting profiles are, of course, different in each case. We permit even point sources and the use 
of thermostats to regulate temperature in different parts of the domain. These and other features 
are illustrated in Examples 1-3 in Sect. 2. We also extend the results in lfT6l in a different direction, 
to "closed" systems, i.e., systems in which the total particle number is conserved and system-bath 
interactions are limited to energy exchange (Sect. 4). 

Our second aim is to understand why memory and finite-size effects are so mild in the models 
tested, or, put differently, why a scheme that ignores these effects can perform so well in profile 
predictions. Here we find the key to lie in cell geometry, referring roughly to the geometric relation 
of cell walls and passageways to the rotating disk. To examine these issues in greater depth, 
we consider a broader class of cell designs. Among the geometries studied here, the one that 
represents the greatest departure from the ideas in IfTBTl (see above) is used in Examples 5 and 
6(A) in Sect. 3.3 (respectively Figs.|9ta) andfTOTa)). In these examples, walls between cells are 
broken down altogether leaving a domain with rotating disks inside. To summarize, we find that 
finite-size and memory effects can be significant or insignificant. By understanding the dynamical 
mechanisms behind them, we learn to identify "good" and "bad" geometries. With regard to 
directional bias in particle trajectories due to memory, an issue also addressed in [15], in Sect. 3.1 
we offer a correction scheme that is proven to be effective for systems with "bad geometry." 

Stochastic idealizations. Alongside the Hamiltonian models, the authors of IfTBI introduced a 
class of Markov jump processes intended to be their stochastic idealizations. Our generalizations 
in Sects. 2 and 4 apply also to these stochastic models. We have run many tests, and have found - 
without exception - that the simulations converge easily to predicted values (with, not surprisingly, 
smaller errors and faster convergence rates than for their Hamiltonian counterparts). While this 
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Figure 1: Illustration of ID chain and baths. 



class of stochastic models, sometimes called "random-halves models" (see Il32l0 . are interesting 
in their own right, we have elected to focus on Hamiltonian models in this paper because they are 
both more complex and more physical. 

Related works. For reviews of some general physical and mathematical issues involved in the 
study of nonequilibrium steady states and transport processes, see, e.g., [|6l [T9ll22ll29ll40ll4Tll44ll . 
Much more is known in the context of stochastic models or models with a significant stochastic 
component. Of the considerable literature on that topic, we cite here a few that are most closely 
related to the present study: (l[lia|5l[ni[Il|251IM[3Ql[3llM[36l. The literature on models 
with purely deterministic internal dynamics is more sparse. In addition to the articles already 
cited, ID|7l|9l[I3[I![I!|2D|23|^ main ones we know of. Finally, aside from 

putting [fT6l in the more general contexts to which it belongs, we have taken some small steps in 
this paper to connect geometry and local dynamics to nonequilibrium properties, a connection that, 
in general, remains to be understood. 



1 Review of models and results in [[I6|| 



Since the present paper is an extension of n\M with a focus on Hamiltonian models, we begin with 
a brief review of the models and relevant findings in the Hamiltonian part of lfT6l . 



1.1 Description of models 

Single cell in isolation. Let r C M 2 be a bounded, simply connected planar region with piecewise- 
C 3 boundary. A disk D is placed in the interior of T ; its center is fixed, and it is allowed to rotate 
freely about its center. The state of the disk is given by (9, u), where 9 is the angle (relative to a 
marked reference point) and to its angular velocity. A number of particles move about in the cell. 
The state of a particle is given by (x, v ), where x G T := T \ D is its position and v E M 2 its 
velocity. Particles fly freely with constant velocity until they hit either dT or dD (they do not 
interact with each other directly). A particle that collides with dT reflects elastically: 

vH4) = v ll (to), v\t+)= ~V ± (t ), 

where t is the collision time, v " is the component of the particle's velocity tangential to dT , and 
v 1 - is the normal component. When a particle hits dD , the particle swaps its tangential velocity 
with the disk's angular velocity: 

v ± (4)= -AO- 
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We assume that the masses of the particles and the moment of inertia of the rotating disk are scaled 
so that the kinetic energy of each particle is \v | 2 and the rotational energy of the disk is co 2 . Energy 
is conserved in the interaction rule above. 

Single cell coupled to two heat baths. Suppose now the cell has two openings, which we identify 
with subsets 7^, 7^ C dT . It is useful (though not necessary) to think of T as having a left-right 
symmetry, and 7^ and 7^ as symmetrically placed on the left and right sides of T respectively. 
These two openings are connected to two heat baths that absorb and emit particles. A particle 
absorbed by one of the baths leaves the system forever. Each bath is characterized by two param- 
eters, a temperature T and an injection rate p. A Poisson clock of rate p is associated with each 
bath. When its clock rings, the bath places a new particle into the system at a (uniformly) random 
position x E 7, where 7 = 7^ for the left bath and 7 = 7r for the right bath. The particle is 
assigned a random (I.I.D.) velocity v sampled from the distribution 

ce~^ \v II sin<z>| dv , c— — ^ , (1) 

'IX 



where (3 = 1/T and p E [0, tt] is the angle between v and 7, measured so that p = ~ points into 
r . The actions of the two baths are independent. 

ID-chain coupled to two heat baths. Consider N copies of the cell described above - call them 
C\, C%, • ■ • , Cn - each with two openings 7l and 7^ symmetrically placed. We think of these cells 
as occupying integer sites 1, 2, • • • ,N, and for each i identify the right opening 7^ (Cj) of the ith 
cell with the left opening 7x(Ci+i) of the (i + l)st cell, so that particles are able to pass from cell 
to cell along the chain. Sites and N + 1 are occupied by heat baths which inject particles into 
(and absorb particles from) cells C\ and Cn respectively. The rules of injection are as above. 

Remarks, (i) We note that except for the chain-bath interaction, which is stochastic in nature, the 
dynamics within the chain are purely deterministic and energy-conserving. It is for this reason that 
we refer to these models as Hamiltonian chains. The local dynamics are, in fact, not symplectic at 
collisions between particles and rotating disks, (ii) Even though particles do not interact directly 
with each other, they do interact quite strongly via the disks, and the overall character of the 
dynamics is very strongly nonlinear. 

Equilibrium distributions. The following result characterizes the Hamiltonian chain in thermal 
equilibrium with two equal heat baths. 

Theorem 1 (Eckmann- Young |fT6l ). Consider an N -chain in contact with two baths, both with 
temperature T and rate p : 

( i) Single cells. For N = 1, the system preserves a probability measure pb T,p characterized by 
the following properties: The number k of particles in the cell is a Poisson random variable 
with mean 

r area(T) p 

N s/t 

where \^\ is the length of an opening. The conditional measure ofp T,p on Qk, the state space 
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for a single cell containing k indistinguishable particles, is given by 

k 

Ck exp ( — (3 [to 2 + \vj | 2 ) j d6 dco ■ dx\ ■ ■ ■ dxt ■ dv\ ■ ■ ■ dvt , 

i=i 

where (3 = 1/T and Ck is the normalizing constant. 

(ii) iV-chains. Let pj' p , i = 1, 2, • • • , N, be N copies of p T,p . An N -chain in contact with two 
equal baths held at (T, p) preserves the product measure ®f = i pf ,p . 

Here are a few different ways to measure temperature for a cell with distribution p T ' p : 

(i) The mean disk energy, E[w 2 ], is |T . 

(ii) The mean energy per particle in the cell, E[|t> | 2 ], is T. 

(iii) The expected energy of a particle exiting the cell is |T; its distribution is given by Eq. (OQ). 

Nonequilibrium steady states. Suppose a chain is connected to two unequal heat baths, i.e., 
(Tl,pl) 7^ (Tr,Pr), where Tl and T R are the temperatures of the left and right heat baths re- 
spectively, and /7l and p R are their respective injection rates. It was found numerically in lfT6l 
(though not proved) that irrespective of initial condition, following a transient the system settles 
down to a nonequilibrium steady state determined solely by the quadruple (T L , p L ; T R , p R ). The 
product nature of the equilibrium distribution in Theorem Q] depends crucially on having detailed 
balance within the system (see lfT6ll ). In nonequilibrium steady states, there is no longer detailed 
balance, and the invariant measure cannot be explicitly computed. In particular, the system has 
spatial correlations, and the steady state distribution is not of product form. 

1.2 Scheme for predicting nonequilibrium profiles proposed in Ifl6ll 

In this subsection we review the scheme proposed in |[T6ll for deducing approximate mean disk 
energies, temperatures, and other profiles in chains that are out of equilibrium. This scheme is 
based on the following ideas: 

Assumption 1. The system is close to having zero directional bias. When a particle leaves a cell, 
let Pl denote the probability that it exits to the left and P R the probability that it exits to the right. 
By zero bias, we mean P L = Pr = ~, regardless of the state of the particle at the time it entered 
the cell or the state of the cell while it is there. 

Assumption 2. The system is close to local equilibrium. Given T L , T R , p L , and p R (with (T L , p L ) ^ 
(T R , p R )) and x £ (0, 1), let Pn,[xN] denote the marginal distribution of the invariant measure for 
an iV-chain at site [xN]. We assume there are numbers T = T(x) and p = p(x) independent of iV 
such that for the iV-chain in question, Pn,[xN] is close to p T{ ~ x )^ x ) where p T,p is as in Sect. 2.1. 

The rationale for the first assumption is that when the length of the opening | — ^ | is small, a 
particle will remain inside each cell it visits for a long time, eventually losing memory of past 
events. The second assumption is based on the idea of local thermal equilibrium (LTE) [|22ll25l , 
which, if true for these models, will imply the convergence of Pn,[xN] to p T ^ x ^p^ as iV — > oo. The 
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idea in |[T6ll is to treat the approximations in Assumptions 1 and 2 as though they were exact, and 
to derive under these simplifying assumptions mean energy and particle density profiles. 

Derivation of profiles assuming zero bias and local equilibrium. All means below are taken at 
steady state. For each cell i, let 

Ji = mean number of particles exiting cell i per unit time, 
Qi = mean total energy exiting cell i per unit time. 

By Ji, we refer here to the total number of particles exiting cell i via either one of its two exits, 
and similarly for Qi. Taking the approximation in Assumption 1 to be exact, we get, from the 
conservation of mass and energy, 

Ji = ^{Ji+i + Ji-i) , Qi = ^(Qi+i + Qi-i) , i = 1,2, • • • ,N , 
with boundary conditions 

3 3 
Jo — 2pL , Jn+i = 2pR , Qo = 2p L • -Tl i Qn+i = 2pR • -T R . 

It follows that the quantities Jj and Qi linearly interpolate between the boundary conditions. As 

N — > oo, J[ x n] and Q\ x n] converge (trivially) to 

J(x) = 2p L (l -x) + 2p R x , 
Q(x) = 3p L T L (l - x) + 3p R T R x 

for all x G (0, 1) . Returning to the iV-chain, the local equilibrium assumption together with 
Theorem Q] then tells us that for cell i, we have Qi = Ji - |Tj, leading to the 

{•-t rpf \ 2 Q(x) p L T L (l - x) + p R T R x 
mean temperature profile T(x) = - ■ ——- = . 

3 J(x) p L {l-x)+p R x 

(We think of the ith cell as occupying location Xi = and the baths as located atx E {0,1}.) 
Similarly, we find the 

mean disk energy profile D(x) := lim Efcuj^i] = -T{x) , 

area(T) J(x) 



mean particle density profile k(x) := lim E[k[ xN ]} = 2y/n ■ 



\l\ 2^T{x)' 

Additional profiles, e.g., total kinetic energy in each cell, are also easily deduced. 

Remarks. Mathematically, that the formulas above give approximate profiles for sufficiently long 
chains of cells with sufficiently small openings will follow from the continuity of the observables 
T, D and k in the double limit of zero-bias and infinite volume, the former to be interpreted as 
|7| 4 0. Such limits are far from trivial. Rigorous results are out of reach at the present time. 
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In |[T6ll , simulations were performed for models in which T is given by the area enclosed by 
four overlapping disks with small openings introduced, and results were reported for ID-chains 
with various parameters. It was found that even for relatively small N, e.g., iV = 30, the predic- 
tions above compare surprisingly well against empirical data. This suggests that for the models 
tested, (i) memory effects leading to directional bias in particle movements are indeed negligible, 
and (ii) local equilibrium is, to a good approximation, achieved at relatively small N, i.e., finite 
size effects are insignificant. 

2 Generalization to larger class of models 

In this section, we generalize the results of lfT6l in a variety of directions, including: 

- Higher dimensions. We focus on 2D models; phenomena in higher dimensions are expected 
to be similar. 

- Lattice types. In two or more dimensions, there are different lattice types. We focus here on 
rectangular and hexagonal geometries. 

- Boundary conditions. Dirichlet, periodic, and reflecting boundary conditions are considered. 
We show that our scheme applies even to point sources. 

- Thermostats. By this we refer to the use of external means to maintain or regulate the tem- 
perature of various parts of the system. 

Instead of giving a formal treatment of each of these generalizations, we will present in Sects. 2. 1- 
2.3 below three examples that systematically illustrate the various effects. 

General setup and issues. The following are features common to all the models to be presented. 
We first discuss them so we can focus on the specific points of interest within each example. 

(1) Basic configurations. We start with a lattice of disjoint fixed (non-rotating) disks in M 2 . Fig. [2] 
shows a rectangular array; a hexagonal array is shown in Fig. [5] Particles move in the space 
between these disks, bouncing off them in specular collisions. The packing is relatively dense, 
with narrow channels between adjacent disks to allow the passage of particles. The area in M 2 
between this array of fixed disks is partitioned into identical cells whose walls are circular arcs 
(pieces of boundaries of the fixed disks) alternating with openings called gaps; see Fig. 2. At the 
center of each cell is a rotating disk, with which particles exchange energy following the rules in 
Sect. 1.1. Disk energy will always refer to the kinetic energy of a rotating disk. 

(2) Mixing rates and geometry. We use R and r to denote respectively the radii of the fixed and 
rotating disks. These two numbers together with the distance between the centers of the fixed disks 
will affect the number of collisions a particle experiences before it finds its way out of a cell. More 
collisions will naturally lead to better equilibration within cells. Too many such collisions, on the 
other hand, will delay mixing on a more global scale and slow down convergence to steady state. 
Information on the expected number of collisions can be deduced from Theorem 1 . For example, 
on each visit to a cell, the number of collisions between a particle and the rotating disk is, on 
average, (2nr) / (total length of all gaps) . 

(3) Continuum limits. Our domains of interest are bounded. Formally, we fix a region Acl 2 with 
piecewise smooth boundary, and consider a sequence of homothetically magnified images A n of 
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Figure 2: Illustration of rectangular lattice system. Here, the system consists of a 4 x 3 array of fixed disks 
of radius R, defining among them a total of 6 cells. At the center of each cell is a rotating disk (shaded with 
centers marked) of radius r. A cell T is identified. 

A. At stage n, the fixed disks in the system are those that meet A n . To obtain continuum limits of 
profiles for disk energy, particle counts, etc., we map A n back to A (so that the profiles for different 
n are all defined on A) and take n — > oo. 

(4) Predicted vs. simulated results. In each of the examples below, we will report the percent- 
age discrepancy between the predicted and simulated values of energy, particle count, etc. Our 
predictions are based on exact versions of the two assumptions in Sect. 11.21 which cannot be true 
for finite chains with finite gap sizes, and the simulated values given are computed empirically for 
finite - and only moderately large - systems. 

Remark on numerical simulations. All numerical results in this paper are computed by direct 
Monte Carlo computations, i.e., expectations are estimated by ergodic averages based on long- 
time simulations. We prepare the system by sampling the equilibrium distribution at a uniform 
temperature T init , where T init is usually chosen to be the maximum of the boundary temperatures, 
and evolve the system for 10 10 events (an event is a single particle-disk collision, bath injection, or 
particle exit). Statistical errors are quantified by a standard batch means method. Typical relative 
errors for the simulations in this section are < 2%. For example, for the results in Fig. [3] below, 
the (absolute) standard error for the D profile is estimated to be < 0.05 and the relative error is 
< 0.4% at the 95% confidence level. 

2.1 Example 1. Rectangular lattice, Dirichlet boundary conditions 

Our first example is a 2D rectangular array on a rectangular-shaped domain with Dirichlet bound- 
ary conditions (to be explained). Fixed disks of radius R 6 (0, 1) are placed at even lattice points 
{0, 2, • ■ ■ , 2M} x {0, 2, • • ■ , 2A^} . We think of R as being close to 1, so that each square config- 
uration of 4 fixed disks bounds an open cell with 4 narrow gaps, with gap size 2(1 — R). In the 
center of each cell, i.e., at odd lattice points 6 {1, • • ■ , 2M — 1} x {1, • • ■ , 2N — 1}, we place a 
small rotating disk of radius r. Heat baths are viewed as occupying the strip of cells just outside 
the system, injecting and removing particles through the associated openings. For simplicity, we 
consider only configurations where all baths along the same edge have the same parameters (T, p). 
For example, if the bath along the left edge has parameters (T L , p L ), then particles with energies 
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distributed according to Eq. (Q~|) with T = T L are injected at rate p L into each cell in the leftmost 
column, through its opening on the left. We refer to the boundary condition described above as a 
"Dirichlet" boundary condition. 

Profile prediction. The scheme reviewed in Sect. 1 1 . 21 is easily generalized to the present setting. 
Formulating Assumption 1 in 2D is straightforward: it is equivalent to requiring that the probability 
of exiting any of the 4 openings to be ~ |. The precise meaning of LTE here is as follows: let 
A = [0, a] x [0,6] C K 2 be the fixed domain on which the continuum limit will be taken, and set 
M = [na] and N = [nb] for positive integers n. Then for every (x, y) 6 (0, a) x (0, b), the marginal 
distribution of the cell at site ([nx], [ny]) tends to p T,p for some T = T(x, y) and p = p(x, y) as 
n — > oo. 

Returning to the finite system described in the first paragraph, let J^i and Qk/ denote the 
particle and energy exit rate from the (fc, £)th cell. The particle balance equation is the discrete 
Laplace equation 

Jk,e = - (Jk+i,e + Jk~i,e + Jk,e+i + Jk,e-i) > 

with boundary conditions J = 4pL on the left edge, J = 4/>r on the right edge, etc. The Qk,£ are 
described by the discrete Laplace equation as well, with boundary conditions Q = Ap^ ■ §7l on 
the left, Q = 4p^ • |Tr on the right, and so on. In the continuum limit, these equations tend to the 
Laplace equation with the given boundary values. Solutions for J and Q are then used to compute 
the mean disk energy and mean particle count through 

D(x,y) — -Tix.y) = -■ — -, 

K ,yj 2 K ,yj 3 J{x,y) ' 



n(x,y) = 2y/n 



4 - tt(R 2 + r 2 



J{x,y) 



2(1 -R) A^/T{x-y) 



These formulas are virtually identical to those for ID-chains given in Sect. 1.2, differing only in 
some constants, and can be deduced the same way. The differences in constants are due to the fact 
that each cell here is in contact with 4 (rather than 2) cells/baths. 

Simulation results. We now compare the predicted profiles to the results of direct simulations. 
Fig.[3]shows some simulation results using the following boundary conditions: 

(3to P ,Ptop) = (30,1) 



(T L ,p L ) = (5,8) 



60 x 41 cells 



(7r,pr) = (50,1) 



(-^bottom j Pbottom ) = (io,i) 

Fixed disks have radius R = 0.95 and rotating disks r = 0.2. These numbers ensure that on 
average, each particle experiences « 18 collisions and hits the rotating disk 3 times on each 
pass through a cell. 
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Fig.[3]shows profiles for the jump rate J, the energy exit rate Q, the mean disk energy D, and the 
particle count k. Boundary values have been filled in for ease of visualization. We have found that 
the discrepancy between simulation results and predicted profiles are greatest in the "L-shaped" 
regions at corners along the boundary of the domain. Away from these regions, the discrepancy is 
< 1.5% . That the discrepancy is larger here is clearly due to discontinuous boundary values and 
finite-size effects. 

Before leaving this example, let us check our intuition against the plotted profiles. Observe 
first that panels (a) and (c), i.e., plots of jump rates J and disk energies D, reflect the boundary 
conditions in a straightforward way. The saddle-shaped profile of Q (panel (b)) is also easy to 
explain: it is because along the 4 edges clockwise from left, Q = 4p x |T has ratios 4:3:5:1. 
A factor that influences particle count is that all else being equal, particles tend to gather in low- 
temperature regions. The reason is that in a steady state, the total particle flux into and out of any 
region must be equal. At the lower temperature end, particles move slowly, hence more of them are 
needed to produce the same flux. For the boundary conditions considered here, the accumulation 
of particles on the left due to colder temperatures is exaggerated by the higher injection rate on the 
same side. This has resulted in the huge spread for k in panel (d), from > 120 on the left to < 5 
on the right. If we change the boundary conditions to pl < Pr, the profile for k would be quite a 
bit flatter. This is both predicted and confirmed in simulations. 

Non-Dirichlet boundary conditions. We have obtained similar simulation results for systems 
with periodic and reflecting boundary conditions in the vertical direction. 

2.2 Example 2. System cooled by thermostats 

The main novel feature of this example is the use of thermostats in the present context^] A second 
feature of note is the use of exit-only boundaries, which for convenience we will refer to as open 
boundaries. 

Thermostats. For our purposes, a thermostat is a disk with the property that whenever a particle 
with velocity v hits it, the tangential component v" of the particle's velocity is instantaneously set 
to a random value u sampled from the distribution 



where T* = is the thermostat temperature. The normal component of v is updated according 
to the usual rule: v L (t + ) = — t> ± (t _ ) . Thus, a particle that hits the thermostated disk sufficiently 
many times would acquire the energy distribution for temperature T*. A cell is thermostated if it 
contains a thermostat in place of a rotating disk. 

We consider in this example a cooling system on a finite cylindrical domain with an open 
boundary at one end. More precisely, the domain is defined by a rectangular array of fixed disks 
with a periodic boundary condition in the vertical direction. The right end of the cylinder is put in 
contact with heat baths, with which it has the same type of interaction as in the previous example. 
The left end of the cylinder has an open boundary, meaning particles can leave but no new particles 

'See, e.g., ||8][T3 USED] QUEUED f° r more examples and thorough discussions of thermostats in nonequilibrium 
statistical mechanics. 
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(a) J(x,y) 



(b)Q(x,y) 




Figure 3: Simulation results for the rectangular lattice with Dirichlet boundary conditions. Bath parameters 
are as given in the text. The lattice contains 60 x 41 rotating disks, which are viewed as located on a regular 
grid within the domain [0,3] x [0,2]. The disk radii are R = 0.95 and r = 0.2. Simulation values are 
plotted over interior grid points, and boundary values are added to ease visualization. For all profiles, the 
discrepancy between computed and predicted profiles are 1-1.2% for (x, y) away from "L"-shaped regions 
near the corners. 
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enter through it. This configuration leads to a particle flux across the system, from right to left. As 
the relatively hot particles move through the system, they are cooled (either directly or indirectly) 
by a row of thermostated cells. For simplicity, we assume the thermostated cells occupy exactly 
one row which extends the full length of the cylinder, and all thermostats in these cells are set to a 
fixed temperature T* considerably lower than that of the heat bath. 

Profile prediction. To obtain predictions for the various profiles, we proceed as follows: 

First, one solves for J as before, the only difference being that the boundary conditions for J 

are periodic in the y direction. The open boundary for x = corresponds to setting J(0, y) = 

for all y. Solving the boundary value problem readily yields J. 

To find Q, the left and right boundaries are treated in a similar fashion as for J. To handle 

thermostats, we assume that the array is mapped to a domain [0, a] x [0, b] with {y = 0} and 

{y = b} identified, and that the thermostats are at y = 0. We then set, for each (x, y) with 

y G {0,b}, Q(x,y) = |J(x, y)T* , where J(x,y) is as in the last paragraph. This is again a 

boundary value problem, which can be solved to give Q. 

Once J and Q have been computed, profiles for D and k are given by the same formulas as in 

the previous example. 

We point out two reasons why one may expect greater discrepancies between predicted and 
simulated results in the present model than in the Dirichlet case: (1) In setting the "boundary 
values" for Q along the row of thermostated cells, we have assumed that the outgoing particles have 
energy |T*. In reality, a highly-energetic particle that enters a thermostated cell is likely to leave 
it with somewhat more energy than |T* : it interacts with the rotating disk only a finite number 
of times, keeping the normal component of its velocity each time. (2) At the open boundary, J 
must necessarily go to zero, as must therefore particle count. Since a healthy number of particles 
is needed in each cell to ensure proper mixing, the situation near this boundary can be anomalous. 

Also, with few particles, a practical concern is inadequate updating during the course of the 
simulation. 

Simulation results. Fig. @] shows simulation results for the following boundary values: 



Open 



Periodic 



Total system size: 60 x 41 cells 



Row 1 disks: Thermstated at T* = 10 



oooooooooooooooooooooo 



(T r ,pr) = (50,10) 



Periodic 



Again, we use R = 0.95 and r = 0.2. Panel (a) shows the mean disk energy D, which is higher on 
the right side due to bath injections and lower along the top and bottom edges due to the cooling 
effects of the thermostats. Panel (b) shows the particle counts. In addition to the cells near the 
baths, particles accumulate near the top and bottom boundaries due to the lower temperatures 
there. Panel (c) compares directly simulation data and predicted values for some cross-sections 
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of the D profile, and panel (d) examines the situation along y = (where the thermostats are) 
for x e [2, 3], i.e., near the bath. As anticipated, the temperature here is higher than its predicted 
values, with the discrepancy rising sharply as x tends to 3. The particle deficit is then consistent 
with the fact that exit rates per particle are higher than predicted. 

The effects highlighted in panel (d), however, seem confined to the thermostated cells. Away 
from them, the discrepancy between predicted and simulation profiles are about 1-4% for the mean 
disk energy and < 7% for the particle count. The latter shows larger discrepancies near the open 
(left) end of the domain, due presumably to n ~ near there. 

2.3 Example 3. Hexagonal array with point sources 

Our third example is a hexagonal array with point sources and open boundaries. The array we 
use consists of fixed disks of radius R placed on the vertices of a triangular mesh; each edge of 
the mesh has length 1 . Each triangular array of 3 fixed disks bounds a cell. At the center of each 
cell, we place a small rotating disk of radius r. See Fig. [5] The system has open boundaries with 
no bath injections. Instead, it is driven by a finite number of point sources which pump particles 
directly into the interior of the array. We say a cell contains a point source with parameters (T#, p*) 
if instead of being an ordinary cell containing a rotating disk, it emits particles having the energy 
distribution in (Q~|) with T = T„ at a rate of p* irrespective of what enters the cell. These particles 
are injected into each of the 3 neighboring cells with equal probability, i.e. at a rate of p*/3 each. 
Notice that by this definition, point sources are "sinks" as well, in that they also absorb particles. 

Profile prediction. Unlike boundary-driven systems, simply letting domain size tend to infinity 
does not lead to a meaningful limit. To see this, consider a single point source of rate p in the center 
of a roughly circular array of cells. In the continuum limit, the jump rate J satisfies the Laplace 
equation with the condition J = on the boundary of the domain and J = p at the location of 
the point source. The only solution is J = almost everywhere. Similar considerations show that 
quantities like mean particle count are also meaningless in this limit. 

Nevertheless, one can still apply the prediction scheme for finite-size systems: we simply solve 
the discrete J and Q balance equations without taking continuum limits @ These equations take the 
form | | 

Jc = q ^ J& and Q c = - ^ Q c , , 

where ~ is the neighbor relation on the hexagonal grid and the equation holds for all cells c that 
do not contain a point source. As before, open boundaries are implemented by imagining a strip of 
cells c just outside the array and setting J c = Q c = for these cells. Finally, a cell c containing a 
point source is implemented by defining the values of J c and Q c for that cell by 

3 

J c = p* and Q c = -p* T* , 

where (T*, p*) are the parameters for the point source in cell c . 

2 Another possibility is to obtain a nontrivial continuum limit by allowing the injection rate p to scale with domain 
size. We do not consider this possibility further. 
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(a) Mean disk energy D(x, y) (b) Mean particle count k(x, y) 




0.5 1 1.5 2 



(c) Slices of D at x = 0.75, 1.5, 2.25 (d) Slices of Q/J and k at y = 



Figure 4: Simulation results for cylindrical array with "cooling" boundary conditions. The array contains 
60 x 41 rotating disks, mapped to grid points inside [0, 3] x [0, 2]. Mean disk energy and particle count 
profiles are shown in (a) and (b). The discrepancy away from edges is about 1-4% in (a) and < 7% in (b). 
In (c), D(x, y) is plotted against y for x = 0.75 (o), x = 1.5 (□), and x = 2.25 (A). The two panels in (d) 
show plots of mean energy per exiting particle Q/J and particle count k along the row of thermostated cells 
at y = 0. Parameters are as described in the text. 
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Figure 5: A patch of hexagonal lattice. 



With the exception of k, which depends explicitly on the geometry of cells, all profiles can be 
derived from J and Q in exactly the same way as before. The expression for k here is easily found 
to be 

v/3/4-7r(iit! 2 + r 2 ) j( xy ) 
2V 71 " ' — 



n(x,y) 



1 - 2R 



Simulation results. Fig. [6] shows the results of a simulation. The domain consists of a hexagonal 
array with an approximately rectangular boundary, with 60 cells along the bottom edge and 40 
along the left. The fixed disks have radii R = 0.485, the rotating disks r = 0.05, so that each 
particle hits the rotating disk on average about 10 times before exiting. We place two point sources 
in the system at the indicated locations. Shown are the disk energy profile D(x, y) and the mean 
particle count k(x, y). The discrepancy between predicted and simulation profiles are < 1.8% for 
D(x, y) and < 3.5% for k(x, y) away from boundaries. 

These results are as expected: the disk energy and particle count profiles clearly reflect param- 
eters of the point sources. One also sees clearly that the number of particles tends to zero as one 
approaches the boundary of the domain, as it should for open boundaries. The sparsity of particles 
explains the relative noisiness near the boundaries of Fig. [(J a) (as discussed in relation to the left 
edge in Example 2). 



Concluding remarks for Section 2 

We have extended the recipe for profile prediction proposed in lfi~6l to a variety of situations and 
tested it in many concrete examples. In each of the three examples presented (and the many others 
not shown), we have found that (i) the predicted profiles are very simple to compute, and (ii) the 
discrepancies between predicted and simulated values are no more than a few percentage points. 
Given that our parameter choices take the systems tested very far out of equilibrium and the system 
sizes considered are only moderate, we find the agreement between predicted and simulated results 
strong. We conclude that the scheme proposed is effective for the class of models studied. 



3 Memory, finite-size effects, and geometry 

The success of the prediction scheme proposed in lfi~6l and extended in Sect. 2 of the present paper 
prompts the following question: Why is it that memory and finite-size effects have so little impact 
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(a) Mean disk energy D(x, y) (b) Mean particle count k(x, y) 



Figure 6: Simulation results for the hexagonal lattice with two point sources. Source A has parameters 
(T, p) = (20, 6); Source B, (T, p) = (100, 70). The domain consists of 60 cells along the top and bottom 
edges and 40 cells along the left and right. (Because of the hexagonal geometry, this means 120 x 40 = 4800 
cells in total.) The discrepancy away from boundaries are < 1.8% for D and < 3.5% for k. 

on macroscopic profiles in the models considered? The purpose of this section is to investigate, in 
the simpler setting of ID chains, which dynamical mechanisms and/or aspects of cell geometry are 
responsible for this outcome, and what happens if the geometry is varied. 

We focus on the following two issues: 

(i) Directional bias due to memory of particle history, i.e., depending on its past history and the 
state of the cell, a particle may have a preferred direction for exiting; and 

(ii) Incomplete equilibration within cells, which may (or may not) happen when particles collide 
too few times with the rotating disk before exiting. 

3.1 Directional bias due to particle history 

The first part of this subsection has some overlap with lfl~5ll . which contains a phenomenological 
study of the same bias questions. Our findings are consistent with those in that earlier study, 
though we take a slightly different approach, and the bias-correction equations we derive (which 
are simpler than those in lfl~5l0 will be tested in models with significant bias. 

The ID-chain used here is exactly one row of the 2D-rectangular configuration used in the 
previous section, with the top and bottom gaps of each cell sealed, i.e., particles reflect off these 
two straight segments elastically. In Parts (A) and (B) below, we consider systems with fixed disk 
radius R = 0.95 (leading to gap size = 0.1) and rotating disk radius r = 0.2, the same parameters 
as those used in Examples 1 and 2 of the previous section. These parameters will be varied, with 
consequences, in Part (C). 

The analysis below applies only to situations when there is a sufficient number of particles in 
the cell. We assume a mean occupation number of at least about 5. 
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(A) Single cell 



We begin by measuring the amount of bias in a single cell using numerical simulations. Consider a 
2x2 array of fixed disks defining a single cell as shown. Heat baths with the following parameters 
are attached to the left and right openings: 



T L = 1 - ±AT l 
p L = 1 - ±Ap" 



Let p _ and p + denote the rates at which particles exit the cell to the left and to the right respectively. 
Similarly, we denote mean energy outflux rates by Q ± , and define the mean "temperatures" of 
outgoing particles by T^ 1 = \Q ± I p ± ■ We express the observed bias by 

AT out = T + _ T - ? Ap out = p+ _ p - _ 

Note that zero directional bias would imply Ap out = . 

We then express AT""' and Ap out in terms of AT" 1 and Ap m , stipulating that T L , T R , and also 
Pl, Pr, be symmetric about 1. Expanding to leading order, we obtain 

AT out w arAT in + p T An™ , 

(2) 

Ap on * « a p AT m + /3 p Ap m , 

where the coefficients are computed (using —0.2 < AT m , Ap m < 0.2) to bed 

a T » 0.130 ±0.001 , /3 T » 0.014 ±0.008 , 
a p « 0.021 ±0.007, P « 0.125 ± 0.001 . 

Relevant dynamical mechanisms. In general, we think of the rotating disk, whose energy fluc- 
tuates wildly, as a randomizing force, causing particles to lose memory of their previous spatial 
locations and directions. Without a doubt, this mechanism is responsible for diminishing bias in 
exit direction^ There is a situation, however, when this mechanism fails: Suppose a particle enters 
from the right and hits the rotating disk at roughly a right angle. If the rotating disk at that moment 
happens to have energy considerably lower than that of the particle, then after the energy exchange 
(see Sect. 1.1), the normal component of the velocity of the particle continues to dominate its 
tangential component, and the particle is reflected back essentially the way it cameJl 

Define a reflected particle to be one which, in the course of its journey through the cell, hits 
only the rotating disk and the two fixed disks on the side of the opening through which it enters, 

3 The error estimates given here are 2x standard error (corresponding to 95% confidence intervals). This degree of 
uncertainty has almost no perceptible effect on the bias correction results in Sect. 3.2(B). 

4 Concave cell walls also help in the scattering of particles, but they appear not to be as important as the randomizing 
effects of the rotating disk. 

5 A similar observation is made in ifTSIl . 




T R = 1 + | AT*' 
Pr = 1 + iAp" 
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i.e. it is "reflected back". (Such particles need not exit the first time around; they can turn around 
and have multiple interactions with the rotating disk.) For AT m = and Ap m G [—0.2, 0.2] , we 
find that among particles entering from the left (respectively, the right), a fraction r ps 0.123 is 
reflected. Since higher energy particles are more likely to reflect, one naturally expects the mean 
temperature of reflected particles to be higher than that of other particles. At (T, p) = (1, 1), we 
find this difference to be t ~ 0.09. 

Assuming that particles that are not reflected are thoroughly randomized, one would expect, in 
terms of orders of magnitude, that a T , f3 p ~ r and f3 T ~ rt. The estimate for /3 T is based on the 
fact that on the side with higher injection rate, a larger fraction of particles is reflected, and this 
fraction has higher temperature. Computed values of cxt , ftp and (3t are consistent with these very 
naive estimates. As for a p , we can only say that it should be > because higher temperature leads 
to more reflected particles. 

The findings above support our contention that reflected particles are an important source of 
bias, but the true story is much more complicated. For example, when a particle interacts with the 
rotating disk, even as some memory is lost, its post-collision trajectory is not entirely random. It 
also retains some memory of its previous energy since the normal component of its velocity upon 
impact is kept. A complete analysis is beyond the scope of this article. 

Linear response for general (T, p). If bias arises mainly from individual particle trajectories — 
and that appears to be the case when there is a sufficient number of particles in the cell — Ap ont 
and AT out should scale linearly with the total injection rate pl + Pr- This allows us to extend the 
linear response relation © to boundary conditions 

(T L , p L ) = (T — i AT m , p - l - Ap m ) , (T R , p R ) = (T + X -AT\ p + \ Ap m ) 
for general (T, p) . The general linear response relations are 

^rpout \rpin Ao^ U 

~' a T ■ — - h Pt 



T T p 

^ out ^ T i n a . 



Up ■ ~^r- + P P , (3) 

p T p 

where the phenomenological coefficients a T , (3 T , etc., are as in Eq. ©. We have numerically 
verified that these scaling relations are accurate to within 2% for 2 < p < 12 and 10 < T < 60@ 



(B) Chains 

A natural next step is to try to use the linear response relation © to account for the effects of 
bias on macroscopic profiles. A priori, knowledge of bias effects in a single cell need not tell us 
about how a cell would behave when placed within a chain: single cells are driven by I.I.D. bath 
injections whereas cells within a chain receive correlated inputs, and long-range correlations are 
known to be substantial in nonequilibrium steady states. 

6 As pointed out in ifTSI . in principle these coefficients may depend on particle density (this follows from dimen- 
sional analysis). We have not found evidence for any significant variation within the range of parameters used here. 
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Prediction scheme with built-in bias correction. We derive phenomenological equations for 
macroscopic profiles assuming that bias at a given site i is determined solely by the temperature, 
particle flux, and their gradients at i, in the same way as in a single cell. Let Qf and Jf be the 
mean rates of energy and particle flux out of cell i, i.e., Qf is the mean total energy that passes 
from cell i to cell i + 1 per unit time, and so on, and let Tf- = \Qf / ■ 
First, we have the balance equations 

pf + PT = pti + pT+i > 
Qf + Q~ = QU + Q m . (4) 

These are augmented by bias corrections of the same form as Eq. ©: 



T+-T7 A7f Ap\ 
- a T ■ — — h Pt ' 



Ti Ti p. 



in 

1 



p+ - P 7 AT in Ap in 

= a p ■ — h p p ■ — , (5) 



T 



where we define 



and 



T i = 2(^-1 + T i+l) ' Pi = ^(Ptl + Pi+l) 



AT™ = T i+1 -T j t 1 , ApT = p- +1 - pti ■ 

Our proposal is to first compute the coefficients a T , ft T , a p , p from single cell simulations, and 
then to solve Eqs. © and (0) numerically to obtain bias-corrected predictions. If a p = f3 p = a T = 
Pt — 0, these equations are precisely the prediction scheme outlined in Sect. 1 1.21 

We have carried out simulations to validate this bias correction scheme. For chains with the cell 
geometry in Part (A), the effect of bias is generally small but measurable. Even with these more 
delicate measurements, our scheme gives quite accurate predictions of Ap out (x) and AT out (x) . 
These findings will be corroborated by those in simulation studies for which bias is more signifi- 
cant, which we now discuss. 



(C) "Good" geometry, "bad" geometry 

Up until now, we have considered a fixed cell geometry (namely that in Fig. [T^a)), which by all 
counts gives rise to very little bias. One would expect bias to increase for geometries which have 
a larger fraction of reflected orbits. For example, one may guess that if one increases the gap size 
I7I = 2(1 — R), it will be easier for reflected particles to exit, and the bias coefficients will go up. 
We have found this to be true, provided I7I is somewhat smaller than 2r. 

A more straightforward way to increase the fraction of reflected particles, however, is to fix gap 
size and increase r, for larger rotating disks are more effective in blocking incoming orbits. Table 
Q] shows the values of cut and Pt for various values of r with gap size fixed at 0.1, equivalently 
with R = 0.95. The trends are self-evident. 

Example 4: A chain with "bad" geometry. We consider here cells with R = 0.95 and r = 0.4; a 
picture of such a cell is shown in Fig.|7£b). This is a candidate for a "bad" geometry: it is almost as 
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r 


0.055 


0.1 


0.2 


0.3 


0.4 


ay 
Pt 


7 x 10" 4 
1.4 x 10~ 3 


0.060 
0.004 


0.13 
0.014 


0.23 
0.023 


0.33 
0.028 



Table 1: Single-cell bias coefficients vs. rotating disk radius. The gap size is \ j\ = 0.1. We observe similar 
trends for a p and j3 p . Below r m j n = 0.05, particles can fly through a cell without being in contact with any 
wall or disk (infinite horizon), and at r max « 0.46, the rotating disk touches the fixed disks. 




(a) R = 0.95 , r = 0.2 (b) R = 0.95 , r = 0.4 (c) R = 0.5 , r = 0.3 

Figure 7: Examples of different cell geometries. Panel (a) shows the geometry used in most examples in 
this paper, including Examples 1-3. The geometry in (b) leads to strong directional bias. In (c), we show a 
geometry that is likely to lead to incomplete equilibration between particles and the rotating disk. 

easy for a particle to bounce back and forth between two adjacent cells as it is to pass from a side 
chamber to the top or bottom chamber of the same cell. Simulations were performed for a chain 
with N = 60, T L = 5, T R = 50, p L = 2, and p R = 12; the results are displayed in Fig. [8] 

Predicted profiles using the scheme without bias correction are shown in dashed lines, while 
bias-corrected profiles are shown as solid lines. As can be seen from the plots, deviations of 
simulation data (circles) from uncorrected profile predictions are nontrivial, unlike the situations 
in Examples 1-3, where cell geometries are considerably more favorable. On the other hand, the fit 
with profiles from our bias -corrected prediction scheme is excellent. These profiles are computed 
from Eqs. © and © using coefficients oct, Pt, ot p , and f3 p from single cell simulations with 
r = 0.4. 

3.2 Incomplete equilibration 

When every particle spends a very long time in a cell before moving on, one may expect the 
particles and the rotating disk in the cell to come close to equilibrating. Thus one may think that 
small gaps (or passageways) between cells speed up the convergence to local equilibrium. In lfT6l 
and [15J, these gaps were taken to be very small: on average, a particle hits the rotating disk more 
than 12 times on each visit to a cell. In Examples 1-3 in the present paper, these conditions are 
relaxed to 3-5 hits per visit. But are small gaps a necessary condition? Put differently, suppose the 
gaps are large (as in Fig. |7jc)), so that each particle moves through every cell quickly, hitting the 
rotating disk only once or twice (or even less) on each pass. Would such incomplete equilibration 
necessarily lead to anomalous behavior? 

We now attempt to address this question. The short answer is: A particle's failure to completely 
equilibrate with its cell environment on each pass is, in and of itself, not an obstruction to normal 
behavior. However, large gaps can amplify finite -size effects. 

For lack of a better term, let us call a cell boundary "porous" if a particle passing through the 



20 







0.2 



0.4 



0.6 



0.8 







0.2 



0.4 



0.6 



0.8 



X 



X 



(a) Mean disk energy D(x) 



(b) Mean particle count k(x) 



Figure 8: Profiles for chains with "bad" geometry. Here, the chains have rotating disk radius r = 0.4 (we 
keep R = 0.95). Shown are simulation data (circles), predicted profile without bias (dashed lines), and 
bias-corrected predictions (solid lines). 

cell is not likely to collide with the rotating disk too many times before exiting. For the type of 
models discussed earlier, porosity is equivalent to larger gap size. We propose below a connection 
between porosity and finite-size effects, and provide numerical confirmation of our thinking. 

Dynamics of equilibration and finite-size effects: a heuristic discussion 

From the point of view of a particle, equilibration takes place as follows: It moves about the chain 
in what is essentially a random walk, exchanging energy with the rotating disks with which it 
collides, hence indirectly with other particles which collide with the same disks. It is an ongoing 
process, a very messy one involving large fluctuations. To say that a particle has "equilibrated" 
basically means that it has acquired "the right statistics." From this picture, it is clear that full 
equilibration in one cell before proceeding to the next is not at all necessary. 

However, a particle has to remain in the system long enough to fully sample the totality of its 
surroundings. For a concrete example, suppose we place a thermostat in the middle of a chain. 
As noted in Example 2, particles retain their normal components in collisions, one might naively 
expect that a particle has to interact with the thermostat many times for equilibration to occur. 
(Actually, it can also acquire this information indirectly.) If the chain is too short and most of the 
particles exit too quickly, then the thermostat may not be fully effective due to insufficient sam- 
pling. In the case of boundary injections, if the domain is not large enough for particles to collide 
a sufficient number of times before exiting the system, then the different sets of statistics carried 
by injected particles (e.g., bath temperature) may not have enough opportunity to get absorbed and 
propagated through the chain. 

The phenomena described above can be thought of as incomplete equilibration with thermostats 
or baths. They will lead to finite-size effects. This kind of equilibration does not require that the 
particle make many collisions in each cell that it visits, though more collisions are conducive to 
better statistics. 

Here is how porosity enters. To exaggerate the problem, consider a situation where cell bound- 
aries are sufficiently porous that on average particles travel ballistically for some (microscopic) 
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distance without exchanging energy with any disks. In such a system, particles will need more 
"room" to make the same number of collisions. The more porous the cell boundaries, the greater 
the mean free path. The discussion above continues to apply, but a longer chain or larger domain 
is required to ensure complete equilibration with thermostats or baths. 

Alternatively, consider a chain of length N in which the mean free path of particles is £. If 
we view particles in this chain as moving on the interval (0, 1), the paths traced out will resemble 
those of a diffusion whose associated lengthscale increases with i (for fixed N) and decreases with 
iV (for fixed £). This suggests the following: (a) Incomplete equilibration in an iV-chain can have 
a smoothing effect on the profiles of relevant observables on (0, 1), as though they are convolved 
with some kernel. In particular, profiles that are strongly curved may lose some of their curvature, 
(b) Whatever the porosity, this smoothing vanishes as iV — > oo. 

The discussion above is entirely heuristic, but we believe it captures the flavor of part of what 
is going on. To summarize, our reasoning implies the following: (1) Other things being equal, the 
more porous the cell geometry, the larger the deviation from expected behavior for a fixed system 
size. (2) In the infinite- volume limit, these effects will vanish, i.e., these are genuinely finite-size 
effects. 

A number of numerical studies were carried out to support (or refute) the thinking above. 
Chains driven by (unequal) boundary injections behave exactly as described: finite-size effects in 
the form of decreased curvature are clearly visible but for the most part not too pronounced for 
chains of moderate length. The thermostat example discussed above provides a more dramatic 
illustration of complete versus incomplete equilibration. 

Example 5: Cells without borders. The basic setup is illustrated in Fig.[9ta): it consists of a ID 
chain connected to two heat baths at the two ends and with a single thermostated cell in (exactly) 
the middle. The boundary conditions are Tl = 80, Tr = 50, and Pl = Pr = 8, and the thermostat 
is set to temperature T* = 10. Neglecting the possible effects of bias, the predictions in Sect. [2]give, 
in the continuum limit, a piecewise linear profile for the disk energy D(x): it linearly interpolates 
the 3 values D(0) = 40, D{\) = 5, and D(l) = 25; particle count k(x) is computed as before, i.e. 
assuming the system is partitioned into cells (with invisible separations), with each cell containing 
exactly one disk. System size will be varied in this study. We consider the following three cell 
geometries: 

- Geometry A (small gaps): R = 0.95, resulting in gap size |7| = 0.1, and r = 0.2; see 
Fig. |7Ja). These are the parameters used in Examples 1-3. The average number of interac- 
tions with the rotating disk is ~ 6 per visit. 

- Geometry B (large gaps): R = 0.5, resulting in gap size I7I = 1, and r = 0.3, as shown in 
Fig. He). This configuration has infinite horizon. Number of hits per pass is « 0.9. 

- Geometry C (no walls between cells): R = and r = 0.5. A chain with this geometry is 
shown in Fig. [9^a). This is as "porous" a geometry as there can be. 

For each of these geometries, we computed D(x) and k(x) for N =31, 61, 121, 201 and 
401 (for Geometry C, note that the derivations in Sect. 11.21 are not affected by the fact that there 
are no cell walls). In each case, deviation from the predicted profile decreases monotonically as 
the length of the chain N increases. Two sets of results for Geometry C are shown in Fig. |9tb). 
As one can see, finite-size effects are very nontrivial at iV = 31, and are still clearly visible at 
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(a) Cells without borders 



N = 31 TV = 401 




(b) Profiles (c) Disk energy distribution, N = 31 



Figure 9: A ID chain without cell borders. A thermostat with temperature T* = 10 is placed at exactly 
x = 0.5; the boundary conditions are Tl = 80, Tr = 50, /?l = PR = 8; r = 0.5. Panel (a) illustrates the 
geometry. In (b), we plot mean disk energy and particle count profiles for two values of N, showing in each 
case both simulation data (curves with open circles) and predictions (solid lines). Note that the value of the 
disk energy profile D(x) at x = 1/2 is fixed by the thermostat to ^T* = 5. Panel (c) shows the energy 
distribution for the 10th rotating disk from the left for N = 31; empirical data (open circles) are plotted 
against the disk energy density for /i T,p with p = 8 and T 44.3, the empirical temperature of the disk. 

iV = 401. From the plotted values of D(x) (set to 5 at x = | on account of the thermostat), one 
sees that the effectiveness of the thermostat is seriously compromised at N = 31: D(x) being 
higher than predicted is consistent with k(x) being lower than predicted. The plot of k shows 
unambiguously the effects of smoothing, confirming the incomplete-equilibration-with-thermostat 
picture proposed. By contrast, the plots for Geometry A (not shown) show no visible finite-size 
effects beyond N = 61; at = 31, the deviations are already smaller than those for Geometry C 
at N = 401. Simulation results for Geometry B are between those of A and C; they contain no 
surprises. 

An intriguing question is the distribution of disk energy for this type of chains, for they are 
the antithesis of chains with tiny gaps (introduced in lfT6ll with the idea that equilibration within 
cells may lead to more rapid convergence to local equilibrium). A specific question is whether as 

— )• oo single site marginals of steady states will tend to measures of the form p T,p defined in 
Sect. 1.1. Numerical evidence suggests that they will: We found that disk energy distributions are 
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ofGibbsformfor even fairly short chains, independently of the completeness of equilibration with 
baths or thermostats. Fig.[9£c) shows the distributions of disk energy at x ~ 0.3 for Geometry C 
with N = 31. The plot is in log-linear scale; data collected from simulations are plotted against 
disk energy distribution for fi T ' p with p = 8 and T equal to the empirical temperature at the site in 
question. Notice the remarkable agreement, demonstrating how close these marginals are to some 
/i T,p long before T approaches its continuum value. 

3.3 Implications for general model designs 

Building on previous ideas, our results in Sects. 2, 3.1, and 3.2 point to a very rich class of out-of- 
equilibrium Hamiltonian models for which we have a simple and effective scheme for predicting 
approximate energy and particle density profiles. We recapitulate here some of the ideas, and 
discuss further implications. 

Retaining the particle-disk interaction introduced in 11331 and ||27ll , we have shown that the ideas 
first put forward in |[T6l are applicable to a very flexible collection of models: 

First, the results are not limited to ID; simulations have been carried out in ID and 2D, and 
we have seen no reason why the type of prediction scheme studied in this paper cannot be applied 
in three and higher dimensions, provided the models have suitable interactions between moving 
particles and pinned-down objects. The type of lattices over which to place the rotating disks 
appears not to be an issue either; we have illustrated it for rectangular and hexagonal lattices. We 
have tested a number of different ways to maintain the system out of equilibrium: be it boundary 
driven (we have tested a variety of boundary conditions), or driven from the interior, using point 
sources or rows of thermostats - in all cases our prediction scheme has produced easy-to-compute 
profiles that show strong agreement with simulations. 

To understand why memory and finite- size effects are so mild in the models tested in IfTBI 
and in Sect. 2, we have proposed an explanation in terms of geometric designs: In these models, 
the cells have mostly concave (scattering) walls, are connected by narrow passageways, and have 
rotating disks at good distances away from the walls. We have shown in Sect. 3.1, as have the 
authors of [fT51 earlier, that such designs give rise to very low degrees of directional bias. We have 
also explained in Sect. 3.2 that the small passageways, which facilitate equilibration within cells, 
also ensure the efficient absorption of bath and thermostat information, reducing finite-size effects. 
This explains to some degree the rather surprising fact that in ifToll . for example, the prediction 
scheme gives remarkably good results in chains with a mere 30 cells. 

Going beyond the models discussed in the last paragraph, we have found bias to be small even 
as we relax considerably the parameters there, i.e., Assumption 1 in Sect. 1.1 applies to a quite 
large collection of models. However, there are equally natural geometries that give rise a non- 
negligible amount of bias. For these geometries, we have derived a correction scheme which is 
relatively straightforward to implement, and shown it to give accurate results (see Sect. 3.1). In a 
different direction, we have produced numerical evidence in support of our ideas on porosity of cell 
boundaries and finite-size effects (Sect. 3.2). If these heuristic ideas could be further developed and 
validated, the implications would be that the results in IfTBI and Sect. 2 of this paper will remain 
valid for systems with large gap sizes (perhaps no walls at all) - provided we go to larger system 
sizes. 

We finish with two examples: Example 6(A), along with Example 5, represents one of the 
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(a) Hexagonal grid without cell borders (b) Mean disc energies vs. x 

Figure 10: Hexagonal model without cell borders. Panel (a) illustrates (to scale) the geometry and placement 
of rotating disks. In (b), we plot mean disk energies, showing simulation data as a function of x (open circles) 
and the predicted profile (solid line); the size of the array is 90 x 3. 

greatest departures from the models studied in |fT6)| for which our prediction scheme remains ef- 
fective; Example 6(B), a close cousin of Example 6(A), makes a connection to the models studied 
in (271. 

Example 6. (A) We consider a cylindrical domain (or rectangular domain with periodic boundaries 
in the vertical direction). Inside this domain, rotating disks are placed on a hexagonal lattice with 
no cell walls or partitions of any kind separating the disks. The disks have radii r = 0.025, and the 
centers of adjacent disks are 1 unit apart. See Fig.fTOla). Heat baths inject particles into the system 
from the left and the right. To deduce energy and particle count profiles, we subdivide the domain 
into hexagonal cells each enclosing exactly one rotating disk in its center, and proceed as before. 
Fig.fTOlb) shows the computed disk energy profile (open circles) and the predicted profile (without 
bias correction) for a system consisting of a 90 x 3 hexagonal array, with (T L , p L ) = (5, 12) and 
(Tr, pr) = (50, 2). A few other sets of parameters besides this one were also tested, with results 
ranging from good to very good. 

(B) The configuration in Fig.fTOla) is reminiscent of the model used in ETl . except that the rotating 
disks in [27| are placed much closer to each other. In (A), we have deliberately kept them apart to 
avoid bias, which for this geometry can be large: We have carried out a simple experiment similar 
to that in Sect. 3.1(A) involving 3 rotating disks of radius 0.475, separated by 3 narrow channels of 
width 0.05. When all 3 disks are thermostated at T = 1, simulations show that a particle entering 
the triangular region between the disks through one of the channels is 3-4 times more likely to exit 
from the same channel than either one of the other two channels. This suggests that a system such 
as the one in (A) but with disks packed this close will exhibit significant bias. Preliminary testing 
confirms this thinking. For example, when p^ = pr, the profile for J, which in the absence of bias 
would be constant along the length of the cylinder, was found to be nontrivially curved (concave). 

Remark. Finally, we observe that as far as methods and overall ideas are concerned, no aspect of 
our studies is tied to the particle-disk interaction used. Thus it is reasonable to conjecture that they 
may extend to other suitable short-range interactions. 
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4 Systems with constant particle number 



Up to this point in the paper, we have studied nonequilibrium steady states in systems where both 
particles and energy can flow across the boundary, with bath injection rates pl and pr playing the 
role of chemical potentials. In this section, we consider systems in which the total particle count is 
maintained at a constant number, and only energy is allowed to flow into and out of the system. The 
difference between these particle-conserving chains and the chains we studied earlier is akin to the 
difference between canonical and grand canonical ensembles in equilibrium statistical mechanics. 

To keep the number of particles fixed in models of the type studied in lfi~6l and in this paper, one 
way is to replace particles upon their reaching the boundary of the domain with new particles whose 
energies reflect those of the baths. That is to say, each particle that exits is replaced instantaneously 
by a new one, and there is no other infusion of particles. Another possibility is to isolate the 
system entirely from its surroundings, and to place in different parts of the domain thermostats set 
to unequal temperatures. Below, we report on findings using particle-conserving baths; we have 
found that simulations using thermostats yield similar results. 

Example 7: ID-chain with particles replaced on exit 

We consider a ID iV-chain as before, with cells i — 1 and i = N being in contact with heat baths. 
Particles are replaced immediately upon exit. To be definite, we assume that when a particle exits 
the system at location x E 7, where 7 is either the left opening of cell 1 or the right opening of 
cell N, independently of its exit velocity it is replaced by a new particle which appears in the same 
position x and has a new velocity v 1 sampled from the bath distribution in Eq. (OQ). Equivalently, one 
can think of the gaps leading outside the chain as being sealed off by two "thermostated walls," 
and when a particle hits one of these walls, it is reflected but with a new velocity. With these 
particle-conserving dynamics, the total number of particles M remains constant in time. Taking a 
continuum limit means letting M, N — > 00 while keeping the density M/N constant. 

Profile prediction. We submit that the various profiles can still be predicted following the basic 
recipe in Sect. \l.2[ one first solves for the mean particle jump rate J(x) and the mean energy 
outflux Q(x); the other profiles are then expressed in terms of J and Q as before. 

What is different here is that because particles are replaced upon exit, the boundary conditions 
for J are Jo = Ji an d Jn = Jn+i ■ These (discrete) homogeneous Neumann boundary conditions 
lead, assuming zero bias, to J(x) = J , i.e., the J profile is constant in space. The Q profile has 
the same boundary conditions as before and is linear. This leads to a linear temperature profile. 
The solutions to the J and Q equations are unique only up to multiplicative constants: they are 
scaled to match the total particle number in the chain. 

Simulation results. Fig.fTTTa) shows simulation results for a chain with N = 60 and M = 2000, 
and boundary conditions Ti = 5 and Tr = 50. Shown are simulation data (dots) and predicted 
profiles - without any built-in bias correction. As can be seen, the agreement is very good except 
for the J profile. 

The discrepancy in the J profile is, at first, a little puzzling: it is somewhat larger than in similar 
non-particle-conserving chains, and is not improved by doubling the total number of particles to 
4000. At the same time, estimates along the lines discussed in Sect. 3.1 show that the amount of 
bias here (Ap out / J < 0.0025) is relatively small. A possible explanation is that in the present 
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(a) Simulation and predicted profiles without bias correction 




X 

(b) Bias-corrected J profile 

Figure 11: Profiles for particle-conserving systems. The system consists of 60 cells and contains 2000 
particles. Shown are simulation data (circles) and predicted profiles (solid lines). The predictions in (a) are 
without bias correction; (b) uses the bias-corrected scheme. Boundary conditions are Tl = 5 and Tr = 50. 



27 



setting, bias can and does lead to a net gradient in J across the length of the chain, whereas in 
non-particle-conserving systems the values of J at the two ends are fixed by the baths. 

Fig. dUb) shows bias-corrected predictions for J using the same numerical values of the bias 
coefficients a p , f3 p , etc., as in Sect. [3l Agreement with simulation data is good, and lends further 
support to the idea that bias is a "local" phenomenon, governed mainly by the mean values of 
thermodynamic quantities and their gradients. 

Conclusions 

The main findings of this paper are: 

(i) For a large class of Hamiltonian systems, macroscopic profiles such as those for energy and 
particle count can be quite accurately predicted in a simple way. The scheme introduced 
in lfl6l and applied to ID chains is extended here to more general settings, including 2D 
systems on rectangular and hexagonal lattices driven by Dirichlet and other boundary condi- 
tions, even point sources. Strong agreement with simulation data is also obtained for systems 
regulated by thermostats, and for systems in which particle numbers are held constant. 

(ii) Memory effects resulting in directional bias of particle trajectories can vary from negligible 
to significant depending on cell geometry. Where such bias is weak, as is the case for a range 
of examples, the prediction scheme above (which assumes zero bias) gives accurate results. 
For "bad" geometries, a scheme with built-in bias correction is shown to be effective. 

(iii) Porous cell boundaries (or large gaps) cause finite-size effects to become more pronounced 
due to incomplete equilibration with baths or thermostats but do not otherwise appear to 
impact macroscopic profiles or LTE. 

A more detailed discussion of (ii) and (iii) is given in Sect. 3.3. 

Acknowledgement. We are grateful to the referee for many detailed and helpful comments. 
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